/* Test Result

Test Size : 500000
Common Product Time:     0.009603
Parallel Product Time:     0.000968
test size : 2000*3000 3000*4000
Matrix Multiplication Time:     5.879109
block mul time :     6.722336
mat = block mat!

*/

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <stdint.h>
#include <string.h>
#include <sys/time.h>
#include <omp.h>

float vout1,vout2;
struct timeval t1, t2;
double timeuse;

// Common Product
void product(int n, float *xs, float *ys)
{
    gettimeofday(&t1, NULL);
    float v0[n];
    for (int i1 = 0; i1 < n; i1 += 1)
    {
        v0[i1] = (xs[i1] * ys[i1]);
    }
    float v2;
    v2 = 0;
    for (int i3 = 0; i3 < n; i3 += 1)
    {
        v2 = (v0[i3] + v2);
    }
    vout1 = v2;
    gettimeofday(&t2, NULL);
    timeuse = (t2.tv_sec - t1.tv_sec) + (double)(t2.tv_usec - t1.tv_usec) / 1000000.0;
    printf("Common Product Time: %15f\n", timeuse);
}

// Parallel Product
void product_parallel(int n,int m,int p, float *xs, float *ys)
{
	gettimeofday(&t1, NULL);
	float v0[n][m];
#pragma omp parallel for
	for (int i1=0; i1 < n; i1 += 1) 
	{ 
#pragma omp parallel for
		for (int i2=0; i2 < m; i2 += 1) 
		{	 
			float v3;
			v3 = 0 ;
			for(int i4=0; i4<p; i4 += 1) 
			{
				v3 = ((xs[i1*m*p+i2*p+i4] * ys[i1*m*p+i2*p+i4]) + v3) ;
			}
			v0[i1][i2] = v3 ;
		}
	}
	float v5;
	v5 = 0 ;
	for(int i6=0; i6<n*m; i6 += 1) 
	{ 
		v5 = (v0[i6%n][i6/n] + v5);
	}
	vout2 = v5;
	gettimeofday(&t2, NULL);
    timeuse = (t2.tv_sec - t1.tv_sec) + (double)(t2.tv_usec - t1.tv_usec) / 1000000.0;
    printf("Parallel Product Time: %15f\n", timeuse);
}

void matrix_mul(int m,int n, int p, float **ma, float **mb, float **out)
{
	gettimeofday(&t1, NULL);
#pragma omp parallel for
	for (int i0=0; i0 < m; i0 += 1) 
	{
		for(int i1=0; i1<p; i1 += 1) 
		{
			float v2[n];
			for(int i3=0; i3<n; i3 += 1) 
			{
				v2[i3] = (ma[i0][i3] * mb[i1][i3]) ;
			}
			float v4;
			v4 = 0 ;
			for(int i5=0; i5<n; i5 += 1) 
			{
				v4 = (v2[i5] + v4) ;
			}
			out[i0][i1] = v4 ;
		}
	}
	gettimeofday(&t2, NULL);
    timeuse = (t2.tv_sec - t1.tv_sec) + (double)(t2.tv_usec - t1.tv_usec) / 1000000.0;
    printf("Matrix Multiplication Time: %12f\n", timeuse);
}

void block_matrix_mul(int m1, int n1, int p1, int m2, int n2, int p2, float ****xs, float ****ys, float ****out1)
{
	gettimeofday(&t1, NULL);
#pragma omp parallel for
	for (int i0=0; i0 < m1; i0 += 1) 
	{	 
#pragma omp parallel for
	for (int i1=0; i1 < p1; i1 += 1) 
	{
		float v2[n1][m2][m2];
		for(int i3=0; i3<n1; i3 += 1) 
		{
			for(int i4=0; i4<m2; i4 += 1) 
			{
				for(int i5=0; i5<m2; i5 += 1) 
				{
					float v6[n2];
					for(int i7=0; i7<n2; i7 += 1) 
					{
						v6[i7] = (xs[i0][i3][i4][i7] * ys[i1][i3][i5][i7]) ;
					}
					float v8;
					v8 = 0 ;
					for(int i9=0; i9<n2; i9 += 1) 
					{
						v8 = (v6[i9] + v8) ;
					}
					v2[i3][i4][i5] = v8 ;
				}
			}
		}
		float v10[m2][m2];
		for(int i11=0; i11<m2; i11 += 1) 
		{
			for(int i12=0; i12<m2; i12 += 1) 
			{
				v10[i11][i12] = 0 ;
			}
		}
		for(int i11=0; i11<n1; i11 += 1) 
		{ 
#pragma omp parallel for
			for (int i12=0; i12 < m2; i12 += 1) 
			{
				for(int i13=0; i13<m2; i13 += 1) 
				{
					v10[i12][i13] = (v2[i11][i12][i13] + v10[i12][i13]) ;
				}
			}
		}
#pragma omp parallel for
		for (int i14=0; i14 < m2; i14 += 1) 
		{
			for(int i15=0; i15<m2; i15 += 1) 
			{
				out1[i0][i1][i14][i15] = v10[i14][i15] ;
			}
		}
	}
	}
	gettimeofday(&t2, NULL);
    timeuse = (t2.tv_sec - t1.tv_sec) + (double)(t2.tv_usec - t1.tv_usec) / 1000000.0;
    printf("block mul time : %12f\n", timeuse);
}


void vector_rand(int n, float *A)
{
    for (int i = 0; i < n; i++)
    {
        A[i] = 1.0 * rand() / RAND_MAX * 100;
    }
}


void matrix_rand(int m, int n, float **A)
{
    for (int i = 0; i < m; i++)
    {
        for (int j = 0; j < n; j++)
        {
            A[i][j] = 1.0 * rand() / RAND_MAX * 100;
        }
    }
}
void mat_block_eq(int m1, int n1, int m2, int n2, float **A, float ****B)
{
    int flag = 1;
    // out =? out1
    for (int i = 0; i < m1; i++)
    {
        for (int j = 0; j < n1; j++)
        {
            for (int k = 0; k < m2; k++)
            {
                for (int t = 0; t < n2; t++)
                {
					float x = B[i][j][k][t] - A[i * m2 + k][j * m2 + t];
                    if ( x < -0.000001 || x > 0.000001)
                        flag = 0;
                }
            }
        }
    }
    if (flag = 1)
    {
        printf("mat = block mat!");
    }
    else
    {
        printf("mat <> block mat!");
    }
}


void mat2block(int m1, int n1, int m2, int n2, float **A, float ****A1)
{
    for (int i = 0; i < m1; i++)
    {
        for (int j = 0; j < n1; j++)
        {
            for (int k = 0; k < m2; k++)
            {
                for (int t = 0; t < n2; t++)
                {
                    A1[i][j][k][t] = A[i * m2 + k][j * m2 + t];
                }
            }
        }
    }
}

int main()
{
    float *va, *vb;
    int l = 500000;
	printf("Test Size : %d\n",l);
	va = (float *)malloc(l * sizeof(float));
    vb = (float *)malloc(l * sizeof(float));
	
    vector_rand (l,va);
	vector_rand (l,vb);
    product(500000,va,vb);
    product_parallel(5,100,1000,va,vb);
	
	int m = 2000;
    int n = 3000;
    int p = 4000;

    float **A,**B,**out,**out1;

    A = (float **)malloc(m * sizeof(float *));
    for (int i = 0; i < m; i++)
    {
        A[i] = (float *)malloc(n * sizeof(float));
    } 

    B = (float **)malloc(p * sizeof(float *));
    for (int i = 0; i < p; i++)
    {
        B[i] = (float *)malloc(n * sizeof(float));
    } 

    out = (float **)malloc(m * sizeof(float *));
    for (int i = 0; i < m; i++)
    {
        out[i] = (float *)malloc(p * sizeof(float));
    } 
	
	printf ("test size : %d*%d %d*%d\n",m,n,n,p);
    

    matrix_rand(m,n,A);
    matrix_rand(p,n,B);
    matrix_mul(m,n,p,A,B,out);

	int m1 = 20;
    int n1 = 30;
    int p1 = 40;
    int m2 = 100;
    int n2 = 100;
    int p2 = 100;
    float ****A1,****B1,****out2;

    A1 = (float ****)malloc(m1 * sizeof(float ***));
    for (int i = 0; i < m1; i++)
    {
        A1[i] = (float ***)malloc(n1 * sizeof(float **));
        for (int j = 0; j < n1; j++)
        {
            A1[i][j] = (float **)malloc(m2 * sizeof(float *));
            for (int k = 0; k < m2; k++)
            {
                A1[i][j][k] = (float *)malloc(n2 * sizeof(float));
            }
        }
    }
	
    B1 = (float ****)malloc(p1 * sizeof(float ***));
    for (int i = 0; i < p1; i++)
    {
        B1[i] = (float ***)malloc(n1 * sizeof(float **));
        for (int j = 0; j < n1; j++)
        {
            B1[i][j] = (float **)malloc(p2 * sizeof(float *));
            for (int k = 0; k < p2; k++)
            {
                B1[i][j][k] = (float *)malloc(n2 * sizeof(float));
            }
        }
    }

    out2 = (float ****)malloc(m1 * sizeof(float ***));
    for (int i = 0; i < m1; i++)
    {
        out2[i] = (float ***)malloc(p1 * sizeof(float **));
        for (int j = 0; j < p1; j++)
        {
            out2[i][j] = (float **)malloc(m2 * sizeof(float *));
            for (int k = 0; k < m2; k++)
            {
                out2[i][j][k] = (float *)malloc(p2 * sizeof(float));
            }
        }
    }



    mat2block(m1,n1,m2,n2,A,A1);
    mat2block(m1,n1,m2,n2,B,B1);
    block_matrix_mul(m1,n1,p1,m2,n2,p2,A1,B1,out2);
    mat_block_eq(m1,p1,m2,p2,out,out2); 
    return 0;
    
}